#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <sys/time.h>
#include <pthread.h>

struct BoxMullerState
{
        double x1, x2, w, y1, y2;
        int useLast;
        struct drand48_data random;
};

void initBoxMullerState(struct BoxMullerState* state)
{
        state->random.__init = 0;
        state->useLast = 0;
        
        struct timeval now;
	gettimeofday(&now, NULL);
	state->random.__x[0] = now.tv_usec;
}

double boxMullerRandom(struct BoxMullerState* state)
{
        double randomNumber;

        if (state->useLast)
        {
                state->y1 = state->y2;
                state->useLast = 0;
        }
        else
        {
                do
                {
                        drand48_r(&state->random, &state->x1);
                        state->x1 = 2.0 * state->x1 - 1.0;
                        drand48_r(&state->random, &state->x2);
                        state->x2 = 2.0 * state->x2 - 1.0;
                        state->w = state->x1 * state->x1 + state->x2 * state->x2;
                }
                while (state->w >= 1.0);

                state->w = sqrt((-2.0 * log(state->w)) / state->w);
                state->y1 = state->x1 * state->w;
                state->y2 = state->x2 * state->w;
                state->useLast = 1;
        }

        return state->y1;
}

//struct que carregara os valores que serao utilizados na paralelizacao
typedef struct black_scholes
{
	double S, E, r, sigma, T; 
	long int m;
	double *trials;
	double means;
}black_scholes;

/*defines da quantidade de threads e struct global
  declaracao tbm da variavel mutex
*/
#define NUM_THREADS 4

black_scholes bs;
pthread_mutex_t mutexbs;


void *calcBlack_Scholes(void *arg)
{
	long int i, start, end, len;
	long offset = (long)arg;
   double random, t;		

	//para criacao dos numeros aleatorios seguindo uma distribuicao normal
	struct BoxMullerState state;	
	initBoxMullerState(&state); 

	//variaveis para controle da area que sera acessado do vetor da struct
	len = bs.m/NUM_THREADS;
	start = offset*len;
	end = start + len;
	
	//variaveis locais para calculo do black_scholes	
	double term1, term2;
	double means = 0, e;                 	
	//algoritmo de black_scholes com monte carlo

	/**
		como cada thread ira acessar posicoes diferentes
		do vetor trials da struct bs, nao eh preciso utilizar o mutex, 
		pois nao havera concorrencia de variaveis
	*/
	for (i = start; i < end; i++)
	{
		random = boxMullerRandom(&state);
		term1 = (bs.r-0.5*pow(bs.sigma, 2))*bs.T;
		term2 = bs.sigma*sqrt(bs.T)*random;
		e = exp(term1 + term2);
		t = bs.S*e;
		if((t-bs.E) <= 0)
			bs.trials[i] = 0;
		else
			bs.trials[i] = exp((-bs.r)*bs.T)*(t-bs.E);
		//calculando a media para aproveitar o laco
		means += bs.trials[i];
	} 	

	pthread_mutex_lock (&mutexbs);
	bs.means += means;
   pthread_mutex_unlock (&mutexbs);	     

	pthread_exit(NULL);	
}

//metodo para calcular o desvio padrao dos dados
double dev(double *v, double mean, long int n)
{
   long int i;
	double sum = 0;
	for(i = 0; i < n-1; i++)
	{
		sum += pow(v[i]-mean, 2);
	}

	return sqrt(sum/n);
}

int main()
{
   pthread_t thread[NUM_THREADS];	
	void *status;
   //variaveis de entrada do modelo de black_scholes  
	scanf("%lf ", &bs.S);
	scanf("%lf ", &bs.E);
	scanf("%lf ", &bs.r);
	scanf("%lf ", &bs.sigma);
	scanf("%lf ", &bs.T);
	scanf("%ld ", &bs.m); 

	//alocando o vetor trials
	bs.trials = (double*)calloc(bs.m, sizeof(double));
	if(bs.trials == NULL)
	{
		printf("ERRO!! FALHA NA ALOCACAO DE MEMORIA\n");
		exit(-1);
	}

	//inicializando o valor da media com zero
	bs.means = 0;
	/*
		iniciando o mutex
		criando threads para calcular os valores das trials e a media
	*/
	pthread_mutex_init(&mutexbs, NULL);
	int i;
	for(i = 0; i < NUM_THREADS; i++)
	{
		pthread_create(&thread[i], NULL, calcBlack_Scholes, (void*)i);	
	}

	//esperando pelas threads
	for(i = 0; i < NUM_THREADS; i++)
	{
		pthread_join(thread[i], &status);
	}


	//calculando o intervalo de confianca
	bs.means = bs.means/bs.m;
	double var = dev(bs.trials, bs.means, bs.m);

	double confwidth = 1.96*var/sqrt(bs.m);
	double confmin = bs.means - confwidth;
	double confmax = bs.means + confwidth;

  //printando os resultados 

	printf("S %lf\n", bs.S);
	printf("E %lf\n", bs.E);
	printf("r %lf\n", bs.r);	
	printf("sigma %lf\n", bs.sigma);
	printf("T %lf\n",bs.T);
	printf("M %ld\n",bs.m);
	printf("Confidence interal: (%lf, %lf)", confmin, confmax);

	//fechando o mutex
   pthread_mutex_destroy(&mutexbs);
	pthread_exit(NULL);
}

